﻿/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_MATH_H
#define _GCTL_MATH_H

// library's head files
#include "../core/enum.h"
#include "../core/macro.h"
#include "../core/array.h"
#include "../core/spmat.h"
#include "math_t.h"

// system's head files
#include "cmath"
#include "random"

namespace gctl
{	
	/**
	 * @brief      符号函数
	 *
	 * @param[in]  a     输入值
	 *
	 * @return     大于零返回1 否则返回-1 等于0则返回0
	 */
	double sign(double a);
	/**
	 * @brief      二分法求一个正数的n次方根
	 *
	 * @param[in]  val    输入值
	 * @param[in]  order  开方次数
	 * @param[in]  eps    终止精度
	 *
	 * @return     开方值
	 */
	double sqrtn(double val, int order, double eps = 1e-5);

	/**
	 * @brief      计算一个椭圆在不同位置的半径
	 *
	 * @param[in]  x_len  x方向半径
	 * @param[in]  y_len  y方向半径
	 * @param[in]  arc    计算方向与x轴正方向的夹角（弧度） 逆时针为正
	 * @param[in]  x_arc  x轴正方向绕原点逆时针旋转的角度 （弧度）
	 *
	 * @return     半径值
	 */
	double ellipse_radius_2d(double x_len, double y_len, double arc, double x_arc = 0.0);

	/**
	 * @brief     计算已椭圆为基准面的高程数据的绝对坐标位置。
	 * 
	 * 假设高程数据的测量方向（如大地水准面）为椭圆切线的垂直方向（与切点与坐标原点的连线方向不一致）
	 * 则高程点的绝对空间位置的维度值与切点的维度值也不一致。此函数可以计算校正后的维度值与球心半径。
	 * 
	 * @param x_len     椭圆的x方向半径（一般为长轴）
	 * @param y_len     椭圆的y方向半径（一般为短轴）
	 * @param arc       计算方向与x轴正方向的夹角（弧度） 逆时针为正（等于纬度）
	 * @param elev      切点的高程值
	 * @param out_arc   校正后的维度值
	 * @param out_rad   校正后高程点的球心半径
	 * @param x_arc     x轴正方向绕原点逆时针旋转的角度 （弧度），默认为0
	 */
	void ellipse_plus_elevation_2d(double x_len, double y_len, double arc, double elev,
		double &out_arc, double &out_rad, double x_arc = 0.0);

	/**
	 * @brief      椭球或者球在不同球面位置的半径
	 *
	 * @param[in]  x_len  x方向半径
	 * @param[in]  y_len  y方向半径
	 * @param[in]  z_len  z方向半径
	 * @param[in]  phi    计算方向与x轴正方向的夹角（弧度） 逆时针为正
	 * @param[in]  theta  计算方向与z轴正方向的夹角（弧度） 逆时针为正
	 *
	 * @return     半径值
	 */
	double ellipsoid_radius(double x_len, double y_len, double z_len, double phi, double theta);

	/**
	 * @brief      一维线性插值
	 *
	 * @param[in]  x1    x1值
	 * @param[in]  x2    x2值
	 * @param[in]  y1    y1值
	 * @param[in]  y2    y2值
	 * @param[in]  x     待插值点坐标
	 *
	 * @return     插值结果
	 */
	double line_interpolate(double x1, double x2, double y1, double y2, double x);

	/**
	 * @brief 一维线性插值 假设数据点是均匀分布的
	 * 
	 * @param xmin x的最小值
	 * @param xmax x的最大值
	 * @param y_arr 一维数据数组
	 * @param x 插值点
	 * @return 插值点的数值
	 */
	double line_interpolate(double xmin, double xmax, const array<double> &y_arr, double x);

	/**
	 * @brief      球面双线性插值函数（单位为弧度）
	 * 
	 * @warning 注意纬度为余纬度
	 * 
	 * lat
	 * |
	 * |
	 * h21----h22
	 * |       |
	 * |       |
	 * h11----h12----> lon
	 *
	 * @param[in]  CoLat1  较大的余纬度
	 * @param[in]  CoLat2  较小的余纬度
	 * @param[in]  Lon1    自西向东较小的经度
	 * @param[in]  Lon2    自西向东较大的经度
	 * @param[in]  CoLat0  待插值点的余纬度
	 * @param[in]  Lon0    待插值点的经度
	 * @param[in]  h11     h11点值
	 * @param[in]  h12     h12点值
	 * @param[in]  h21     h21点值
	 * @param[in]  h22     h22点值
	 *
	 * @return     插值点值
	 */
	double sph_linear_interpolate(double CoLat1, double CoLat2, double Lon1, double Lon2, 
		double CoLat0, double Lon0, double h11, double h12, double h21, double h22);

	/**
	 * @brief      球面双线性插值函数 以度为单位的版本
	 * 
	 * @warning 注意纬度为余纬度
	 * 
	 * lat
	 * |
	 * |
	 * h21----h22
	 * |       |
	 * |       |
	 * h11----h12----> lon
	 *
	 * @param[in]  CoLat1  较大的余纬度
	 * @param[in]  CoLat2  较小的余纬度
	 * @param[in]  Lon1    自西向东较小的经度
	 * @param[in]  Lon2    自西向东较大的经度
	 * @param[in]  CoLat0  待插值点的余纬度
	 * @param[in]  Lon0    待插值点的经度
	 * @param[in]  h11     h11点值
	 * @param[in]  h12     h12点值
	 * @param[in]  h21     h21点值
	 * @param[in]  h22     h22点值
	 *
	 * @return     插值点值
	 */
	double sph_linear_interpolate_deg(double CoLat1, double CoLat2, double Lon1, double Lon2, 
		double CoLat0, double Lon0, double h11, double h12, double h21, double h22);

	/**
	 * @brief      规则网络插值 长方形内数据插值
	 * 
	 * @note 距离平方反比
	 * 
	 * y
	 * |
	 * |
	 * d3-----------d2
	 * |            |
	 * |            |
	 * |            |
	 * d0-----------d1--->x
	 * 左下角坐标x0 y0
	 * 块体尺寸dx dy
	 * 插值点坐标x y
	 *
	 * @param[in]  x0    左下角坐标x0
	 * @param[in]  y0    左下角坐标y0
	 * @param[in]  dx    块体尺寸dx
	 * @param[in]  dy    块体尺寸dy
	 * @param[in]  x     待插值点x
	 * @param[in]  y     待插值点y
	 * @param[in]  d0    d0值
	 * @param[in]  d1    d1值
	 * @param[in]  d2    d2值
	 * @param[in]  d3    d3值
	 *
	 * @return     插值点值
	 */
	double rect_interpolate(double x0, double y0, double dx, double dy, double x, double y,
		double d0, double d1, double d2, double d3);

	/**
	 * @brief      规则网络插值 长方体内数据插值
	 * 
	 * @note 距离平方反比
	 * 
	 *        y
	 *       /
	 *      /
	 *     3------------2
	 *    /|           /|
	 *   / |          / |
	 *  0------------1------> x
	 *  |  7         |  6
	 *  | /          | /
	 *  |/           |/
	 *  4------------5
	 *  |
	 *  |
	 *  z
	 *  左上角坐标x0 y0 z0
	 *  块体尺寸dx dy dz
	 *  插值点坐标x y z
	 *  八个角点值
	 *
	 * @param[in]  x0    左上角坐标x0
	 * @param[in]  y0    左上角坐标y0
	 * @param[in]  z0    左上角坐标z0
	 * @param[in]  dx    块体尺寸dx
	 * @param[in]  dy    块体尺寸dy
	 * @param[in]  dz    块体尺寸dz
	 * @param[in]  x     插值点坐标x
	 * @param[in]  y     插值点坐标y
	 * @param[in]  z     插值点坐标z
	 * @param[in]  d0    角点值d0
	 * @param[in]  d1    角点值d1
	 * @param[in]  d2    角点值d2
	 * @param[in]  d3    角点值d3
	 * @param[in]  d4    角点值d4
	 * @param[in]  d5    角点值d5
	 * @param[in]  d6    角点值d6
	 * @param[in]  d7    角点值d7
	 *
	 * @return     插值点值
	 */
	double cube_interpolate(double x0, double y0, double z0, double dx, double dy, double dz, 
		double x, double y, double z, double d0, double d1, double d2, double d3, double d4, 
		double d5, double d6, double d7);

	/**
	 * @brief      球坐标下规则网络插值 Tesseroid 内数据插值
	 * 
	 * @note 首先是线性插值 然后是球面线性插值
	 * 
	 * 六面体示意图
	 *        lat
	 *       /
	 *      /
	 *     7------------6
	 *    /|           /|
	 *   / |          / |
	 *  4------------5------> lon
	 *  |  3         |  2
	 *  | /          | /
	 *  |/           |/
	 *  0------------1
	 *  |
	 *  |
	 *  rad
	 *  左下角坐标lon0 lat0 rad0
	 *  块体尺寸dlon dlat drad
	 *  插值点坐标lon lat rad
	 *  八个角点值
	 *
	 * @param[in]  lon0  左下角坐标rad0
	 * @param[in]  lat0  左下角坐标lat0
	 * @param[in]  rad0  左下角坐标rad0
	 * @param[in]  dlon  块体尺寸dlon
	 * @param[in]  dlat  块体尺寸dlat
	 * @param[in]  drad  块体尺寸drad
	 * @param[in]  lon   插值点坐标lon
	 * @param[in]  lat   插值点坐标lat
	 * @param[in]  rad   插值点坐标rad
	 * @param[in]  d0    角点值d0
	 * @param[in]  d1    角点值d1
	 * @param[in]  d2    角点值d2
	 * @param[in]  d3    角点值d3
	 * @param[in]  d4    角点值d4
	 * @param[in]  d5    角点值d5
	 * @param[in]  d6    角点值d6
	 * @param[in]  d7    角点值d7
	 *
	 * @return     插值点值
	 */
	double cube_interpolate_sph(double lon0, double lat0, double rad0, double dlon, double dlat, 
		double drad, double lon, double lat, double rad, double d0, double d1, double d2, 
		double d3, double d4, double d5, double d6, double d7);

	/**
	 * @brief      按距离反比加权计算均值
	 *
	 * @param      dis_vec  距离向量
	 * @param      val_vec  数值向量
	 * @param[in]  order    距离加权的阶次 默认为1
	 *
	 * @return     加权平均值
	 */
	double dist_inverse_weight(std::vector<double> *dis_vec, std::vector<double> *val_vec, int order = 1);

	/**
	 * @brief      使用cosine函数对一维数组进行外插
	 *
	 * @param      in_arr     输入数组
	 * @param[in]  in_num     数组长度
	 * @param      out_arr    输出数组
	 * @param[in]  out_num    数组长度
	 * @param[in]  end_value  外插后的端点值，若为NAN则函数会计算原数组端点值的均值作为新的端点值
	 */
	void cosine_extrapolate_1d(double *in_arr, int in_num, double *out_arr, int out_num, double end_value = NAN);

	/**
	 * @brief      使用cosine函数对一维数组进行外插
	 *
	 * @param[in]  in_arr   输入数组
	 * @param      out_arr  输出数组
	 * @param[in]  end_val  外插后的端点值，若输入值为NAN函数会计算原数组端点值的均值作为新的端点值
	 * @param[in]  out_len  外插后的数组长度，若输入值小于0函数会计算最小的2的次方值作为输出的数组长度
	 */
	void cosine_extrapolate_1d(const array<double> &in_arr, array<double> &out_arr, double end_val = NAN, int out_len = -1);

	/**
	 * @brief      使用cosine函数对二维数组进行外插
	 *
	 * @param[in]  in_arr    输入数组
	 * @param      out_arr   输出数组
	 * @param      ori_row   原数据在扩边后的起始行号
	 * @param      ori_col   原数据在扩边后的起始列号
	 * @param[in]  end_valx  x方向外插后的端点值，若输入值为NAN函数会计算原数组端点值的均值作为新的端点值
	 * @param[in]  end_valy  y方向外插后的端点值，若输入值为NAN函数会计算原数组端点值的均值作为新的端点值
	 * @param[in]  out_lenx  x方向外插后的数组长度，若输入值小于0函数会计算最小的2的次方值作为输出的数组长度
	 * @param[in]  out_leny  y方向外插后的数组长度，若输入值小于0函数会计算最小的2的次方值作为输出的数组长度
	 */
	void cosine_extrapolate_2d(const _2d_matrix &in_arr, _2d_matrix &out_arr, 
		int &ori_row, int &ori_col, double end_valx = NAN, double end_valy = NAN, int out_lenx = -1, int out_leny = -1);

	/**
	 * @brief      使用cosine函数对二维数组进行外插
	 * 
	 * 按一维数组排列（列优先，按行排列）
	 *
	 * @param[in]  in_arr    输入数组
	 * @param      out_arr   输出数组
	 * @param[in]  in_row    输入数组的行数
	 * @param[in]  in_col    输入数组的列数
	 * @param      out_row   输出数组的行数（如果输入值为负则函数会计算最小的2的次方值作为输出的行数，否则即按照输入值进行扩边）
	 * @param      out_col   输出数组的列数（如果输入值为负则函数会计算最小的2的次方值作为输出的行数，否则即按照输入值进行扩边）
	 * @param      ori_row   原数据在扩边后的起始行号
	 * @param      ori_col   原数据在扩边后的起始列号
	 * @param[in]  end_valx  x方向外插后的端点值，若输入值为NAN函数会计算原数组端点值的均值作为新的端点值
	 * @param[in]  end_valy  y方向外插后的端点值，若输入值为NAN函数会计算原数组端点值的均值作为新的端点值
	 */
	void cosine_extrapolate_2d(const array<double> &in_arr, array<double> &out_arr, int in_row, int in_col, 
		int &out_row, int &out_col, int &ori_row, int &ori_col, double end_valx = NAN, double end_valy = NAN);

	/**
	 * @brief      使用零值对二维数组进行外插
	 * 
	 * 按一维数组排列（列优先，按行排列）
	 *
	 * @param[in]  in_arr    输入数组
	 * @param      out_arr   输出数组
	 * @param[in]  in_row    输入数组的行数
	 * @param[in]  in_col    输入数组的列数
	 * @param      out_row   输出数组的行数（如果输入值为负则函数会计算最小的2的次方值作为输出的行数，否则即按照输入值进行扩边）
	 * @param      out_col   输出数组的列数（如果输入值为负则函数会计算最小的2的次方值作为输出的行数，否则即按照输入值进行扩边）
	 * @param      ori_row   原数据在扩边后的起始行号
	 * @param      ori_col   原数据在扩边后的起始列号
	 */
	void zeros_extrapolate_2d(const array<double> &in_arr, array<double> &out_arr, int in_row, int in_col, 
		int &out_row, int &out_col, int &ori_row, int &ori_col);

	/**
	 * @brief      查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
	 *
	 * @param[in]  in_array    输入数组
	 * @param[in]  array_size  数组大小
	 * @param[in]  in_val      查找值
	 * @param      index       返回的索引值
	 *
	 * @return     成功0失败-1
	 */
	int find_index(const double *in_array, int array_size, double in_val, int &index);

	/**
	 * @brief      查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
	 *
	 * @param      in_array  输入数组
	 * @param[in]  in_val    查找值
	 * @param      index     返回的索引值
	 *
	 * @return     成功0失败-1
	 */
	int find_index(array<double> *in_array, double in_val, int &index);

	/**
	 * @brief      计算一维分形模型
	 *
	 * @param      out_arr     输出数组
	 * @param[in]  l_val       分形计算的左端点值
	 * @param[in]  r_val       分形计算的右端点值
	 * @param[in]  maxi_range  最大变化值
	 * @param[in]  smoothness  变化光滑度
	 */
	void fractal_model_1d(array<double> &out_arr, int out_size, double l_val, 
		double r_val, double maxi_range, double smoothness);

	/**
	 * @brief      计算二维分形模型
	 *
	 * @param      out_arr     输出数组
	 * @param[in]  dl_val      分形计算的左下角端点值
	 * @param[in]  dr_val      分形计算的右下角端点值
	 * @param[in]  ul_val      分形计算的左上角端点值
	 * @param[in]  ur_val      分形计算的右上角端点值
	 * @param[in]  maxi_range  最大变化值
	 * @param[in]  smoothness  变化光滑度
	 */
	void fractal_model_2d(_2d_matrix &out_arr, int r_size, int c_size, double dl_val, 
		double dr_val, double ul_val, double ur_val, double maxi_range, double smoothness);

	/**
	 * @brief      使用牛顿迭代法计算一个矩阵的近似逆
	 * 
	 * 此方法适合用于对角占优的方阵求逆。
	 *
	 * @param[in]  in_mat       输入矩阵
	 * @param      inverse_mat  逆矩阵
	 * @param[in]  epsilon      迭代的终止精度
	 * @param[in]  iter_times   迭代的次数。默认为10
	 * @param[in]  initiate     对逆矩阵进行初始化
	 * 
	 * @return     最终的迭代精度
	 */
	double newton_inverse(const _2d_matrix &in_mat, _2d_matrix &inverse_mat, 
		double epsilon = 1e-8, int iter_times = 10, bool initiate = true);

	/**
	 * @brief      使用牛顿迭代法计算一个矩阵的近似逆
	 * 
	 * 此方法适合用于对角占优的方阵求逆。
	 *
	 * @param[in]  in_mat       输入矩阵（稀疏的）
	 * @param      inverse_mat  逆矩阵
	 * @param[in]  epsilon      迭代的终止精度
	 * @param[in]  iter_times   迭代的次数。默认为10
	 * @param[in]  initiate     对逆矩阵进行初始化
	 * 
	 * @return     最终的迭代精度
	 */
	double newton_inverse(const spmat<double> &in_mat, _2d_matrix &inverse_mat, 
		double epsilon = 1e-8, int iter_times = 10, bool initiate = true);

	/**
	 * @brief      使用施密特正交化将线性无关的向量组a转换为标准正交向量组e
	 * 
	 * @note       经过正交变化的向量组中任意两个向量正交。
	 *
	 * @param[in]  a     个数为a_s的线性无关的向量组的指针（这里将向量组保存为一个连续的一维数组）
	 * @param      e     输出的转换后的标准正交向量组（保存为一个连续的一维数组）
	 * @param[in]  a_s   向量的个数
	 */
	void schmidt_orthogonal(const array<double> &a, array<double> &e, int a_s);

	/**
	 * @brief      一维数组差分（使用二阶差分公式）
	 * 
	 * 计算一个一维数组中相邻元素间的差分结果（求导）。
	 *
	 * @param[in]  in       输入数组
	 * @param      diff     输出的差分结果
	 * @param[in]  spacing  相邻元素的距离
	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。两边的数据将分别使用对应的向前或向后差分公式
	 */
	void difference_1d(const array<double> &in, array<double> &diff, double spacing, int order = 1);

	/**
	 * @brief      二维数组差分（使用二阶差分公式）
	 * 
	 * 计算一个二维数组中相邻元素间的差分结果（求导）。
	 *
	 * @param[in]  in       输入数组
	 * @param      diff     输出的差分结果
	 * @param[in]  spacing  相邻元素对应方向的距离
	 * @param[in]  d_type   求导的类型
	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
	 */
	void difference_2d(const _2d_matrix &in, _2d_matrix &diff, double spacing, gradient_type_e d_type, int order = 1);

	/**
	 * @brief      二维数组差分（使用二阶差分公式）
	 * 
	 * 计算一个二维数组中相邻元素间的差分结果（求导）。数组以列优先的方式储存为一个一维数组
	 *
	 * @param[in]  in       输入数组
	 * @param      diff     输出的差分结果
	 * @param[in]  row_size 数组二维排列的行数
	 * @param[in]  col_size 数组二维排列的列数
	 * @param[in]  spacing  相邻元素对应方向的距离
	 * @param[in]  d_type   求导的类型
	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
	 */
	void difference_2d(const array<double> &in, array<double> &diff, int row_size, int col_size, 
		double spacing, gradient_type_e d_type, int order = 1);
}

#endif // _GCTL_MATH_H